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The  areas  of  research  studied  under  this  contract 
v/ere:  a)  numerical  naethods  for  the  integration  of  stiff 
systems  of  ordinary  differential  equa.tions;  b)  computational 
techniques  for  solving  sparse  systems  of  linear  equations; 
c)  discrete  methods  for  solving  optimal  control  problems 
with  state  space  constraints.  In  addition,  some  prelinjinary 
studies  were  made  to  see  how  the  results  of  the  first  two 
areas  could  be  most  effectively  employed  in  a  design  opti¬ 
mization  pcickage  for  electrical  networks.  At  the  end  are 
a  list  of  papers  tha.t  have  been  published  since  January  1968 
in  outside  journals,  and  another  list  of  papers  with  the 
current  status  on  their  publication  in  a  journal. 
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Numerical  Methods  for  Stiff  Syatems 

The  design  of  efficient  numerical  integration  methods 
for  stiff  systems  of  ordinary  differential  equations  has  been 
studied  further.  The  formulae  designed  earlier  were  all  of 
a  one-step  type.  In  order  to  combine  in  those  one-step 
methods  a  reasonably  large  order  of  accuracy  in  the  conven¬ 
tional  sense  v/ith  strong,  global  stability  and  accuracy  to 
cope  v/itli  stiffness,  one  has  to  consider  formulae  containing 
higher  (in  particular,  second)  derivatives.  One-step 
methods  have  the  advantage  of  not  introducing  parasitic 
solutions,  but  to  express  the  higher  derivatives  on  must 
differentiate  the  differential  equations;  i.  e.  ,  compute  the 
Jacobian  matrix.  Recently,  more  attention  has  been 
focused  on  linear  multistcp  methods  which  provide  for  con¬ 
siderable  flexibility  without  containing  higher  derivatives. 
However,  they  do  introduce  parasitic  roots.  An  order  of 
accuracy  p  >  2  is  not  compatible  with  A -stability  in  the 
Dahlquist  sense.  Thus,  slightly  weaker  stability  concepts, 
which  we  refer  to  as  quasi-A-stability,  must  be  adopted  in 
studying  multistep  formulae;  e.  g.  ,  A (of) -stability  as  defined 
by  Widlund.  A  theorem  has  been  proved  which  leads  to  a 
generalisation  of  the  easy-to-apply  test  criterion  for  A-stability 
to  a  general  class  of  quasi-A-stability  concepts.  The  analysis 
has  been  carried  out  in  detail  for  circular,  angular,  and  para¬ 
bolic  domains  in  the  complex  plane  and  applied  in  particular 
to  the  class  of  backward  differentiation  formulae  used  extens¬ 
ively  by  Gear. 

It  has  also  been  possible  to  generalize  the  concept  of 
exponential  fitting,  introduced  earlier  for  one-step  formulae, 
to  multistep  methods.  This  leads  to  the  definition  of  a  one- 
parameter  family  of  multistep  methods  with  global  accuracy 
properties.  As  special  cases,  this  family  contains  the  well 
known  closed  Adams  formulae  (fitted  at  q  =  0  ),  and  the  back¬ 
ward  differentiation  formulae  (fitted  at  infinity).  It  has  been 
shown,  for  the  class  of  all  two-step  formulae  of  order  p  ^  2 
whose  A-stability  was  stvidied  earlier,  that  exponential  fitting 
is  compatible  with  both  A-stability  and  boundedness  of  the 
local  truncation  error  uniformly  with  respect  to  the  location 
of  fitting. 


Exponential  fitting  ties  a  given  formulae  to  a  particular 
problem  one  is  solving  or,  more  specifically,  to  the  spectrum 
of  the  accompanying  Jacobian  matrix.  In  general,  the  location 
of  fitting  and  the  parameter  in  the  formula  then  vary  with  the 
independent  variable.  As  an  alternative  to  exponential  fitting, 
one  can  attempt  to  make  a  good  a  priori  choice  of  the  para¬ 
meters  contained  in  the  formulae  which  are  then  kept  fixed. 

It  was  shown  that  for  the  one-step,  one -parameter  formula 
considered  in  RC  1852,  there  exists  an  optimal  para¬ 
meter  choice  in  a  certain  sense. 


The  computation?!  methods  for  solving  sparse 
systems  of  linear  equations  were  developed  further.  As 
more  applications  were  examined,  the  need  for  various 
adaptations  or  even  new  methods  was  recognized.  Many 
of  these  were  implemented,  and  are  summarized  below. 

In  addition,  it  was  felt  that  there  was  a  diverse  collection 
of  scientists  working  on  sparse  matrices  who,  because 
they  were  in  different  scientific  disciplines,  did  not  kirow 
what  had  been  develoj^ed  in  other  areas.  A  symposium  on 
sparse  matrix  methods  and  applications  was  organized 
and  held  September  9-10,  1968.  The  Proceedings,  con- 
taijiing  extended  abstracts  of  the  talks  with  complete 
bibliographies  and  the  text  of  a  panel  discussion,  is 
forthcoming,  and  should  serve  to  unify  the  references 
in  this  field, 

A  major  study  was  made  of  direct  methods  for 
solving  linear  equations  as  they  apply  to  sparse  matrices. 
The  two  leading  methods,  judged  by  the  applications, 
were  the  various  forms  of  Gaussian  elimination  and 
Gaus  s -Jordan  elimination.  A  careful  comparison  was 
made  between  the  elimination  form  of  the  inverse  (EFI) 
and  the  product  form  of  the  inverse  (PFI),  respectively 
a  form  of  Gaussian  and  Gauss -Jordan  elimination.  The 
PFI  is  predominantly  used  in  the  major  linear  program¬ 
ming  codes  in  use  today.  One  of  the  results  of  this  study 
is  that  the  EFI  method  is  superior  for  sparse  matrices 
in  that  it  requires  less  storage  and  arithmetic  operations 
than  the  PFI.  The  conclusion  is  that  any  newly  developed 
linear  programming  codes  should  incorporate  the  EFJ. 
method.  To  give  quantitative  substance  to  this  conclusion, 
some  statistics  were  compiled  for  the  computation  time 
and  storage  required  for  both  the  PFI  and  EFI  in  solving 
systems  of  eqviations  where  the  coefficient  matrix  was  a 
sparse  matrix  with  a  random  sparseness  structure. 

A  more  theoretical  investigation  was  made  to 
understand  the  fill-in  phenomenon  and  to  characterize 


algorithms  which  require  the  least  amount  of  storage. 

In  the  paper  this  is  called  s-minimal.  The  fill-in 
phenomenon  refers  to  the  fact  that  when  one  decomposes 
a  sparse  matrix  into  a  form  of  its  inverse,  it  always  re¬ 
quires  more  storage;  i.  e.  ,  the  matrix  fills  in.  It  was 
discovered  that  both  the  EPT  and  PFI  as  well  as  the 
Householder  algorithm  for  factoring  A  into  QU  arc  not 
s-minimal.  This  means  that  certain  computations  and 
storage  locations  apparently  necessary  for  these  algorithms 
are  really  unnecessary  and,  therefore,  waste  storage  and 
computation  time.  However,  it  was  shown  that  if  the 
original  matrix  had  non7.ero  diagonal  elements,  both  the 
EFI  and  the  PFI  are  s-minimal.  This  result  may  have 
implications  for  rearrangement  algorithms  and  explains 
the  occurrence  of  almost  zero  elements  in  the  linear  pro¬ 
gramming  codes. 

Another  area  in  sparse  matrix  theory  that  was 
studied  was  the  different  methods  of  modifying  the  form 
of  the  inverse  when  the  original  matrix  A  is  changed. 

Five  methods,  some  of  them  new,  were  discussed  and 
compared.  One  of  the  methods,  Kron’s  method,  which 
has  been  discussed  extensively  in  electrical  engineering 
literature,  was  related  to  the  method  used  in  linear  pro¬ 
gramming  codes  with  the  PFI,  This  analysis  puts  Kron's 
method  in  a  new  light  and  it  seems  fair  to  say  that  with 
sparse  matrix  codes  available,  Kron's  method  has  little 
to  offer. 

The  status  of  the  various  sparse  matrix  codes'  is 
summarized  below. 

Interpreter.  This  program  is  a  modification  of 
GNSO  (an  algorithm  which  produced  FORTRAN  statements 
for  solving  a  given  sparse  system  of  equations).  Instead  of 
producing  FORTRAN,  it  produces  a  sequence  of  integers 
which  can  be  interpreted  rapidly  to  produce  the  correct 
solution  of  a  given  sparse  system  of  equations.  For  a 
problem  of  zie  199  ^  199  which  was  suggested  previously 
by  Air  Force  personnel,  68,250  integers  were  generated 
in  15  seconds  on  a  360/67.  This  sequence  of  integers 


was  made  resident  on  a  2311  disk  pack  and  had  to  be 
brought  into  core  to  solve  the  problem,  making  the  exe¬ 
cution  I/O  bound.  Even  so,  it  took  only  3  seconds  to 
bring  in  this  code  and  solve  the  199  X  199  system.  The 
CPU  was  idle  1/6  of  the  time. 

Compiler.  A  special  purpose  compiler  was  designed 
and  implemented  to  work  with  GNSO.  This  compiler  trans¬ 
lates  code  gcJierated  by  GNSO  into  IBM  3  60  machine 
executable  code.  The  need  for  such  a  compiler  was  that 
either  the  standard  FORTRAN  compiler  was  too  slow, 
or  could  not  compile  the  size  of  program  generated  by 
GNSO.  For  e:<am23le,  on  a  problem  of  size  57  X  57, 
FORTRAN  H  failed  to  compile  the  FORTRAN  statements 
px'oduced  by  GNSO  because  there  were  too  many  state¬ 
ments,  FORTRAN  G  took  9  minutes,  and  the  new 
sjDccial  compiler  took  1  second.  GNSO  and  the  compiler 
took  29  seconds  to  produce  the  code  for  the  199  X  199 
problem.  This  code  executes  in  .2  seconds. 

New  GNSO.  A  new  version  of  GNSO  based  on  a 
monotone  sequence  of  integers  representing  the  nonzero 
entries  of  a  sparse  matrix  was  implemented.  This  replaces 
the  use  of  the  bit  map  in  GNSO  and  is  more  efficient  for 
large  matrices. 

Band  Matrix  Program.  This  is  a  special  program 
for  sparse  matrices  which  have  a  band  structure;  i.  e. ,  all 
the  nonzero  elements  lie  in  a  band  about  the.  diagonal.  The 
program  is  coded  in  FORTRAN,  but  in  the  parts  of  the 
code  which  are  tho  most  repetitive,  machine  language  is 
used  to  achieve  maximum  efficiency.  For  computers  like 
the  IBM  3  60  models  91  and  95,  the  program  takes  full 
advantage  of  the  loop  mode  so  that  maximum  vitilization  of 
the  CPU  can  be  approached. 

Sparse  Matrix  Gaussian  Elimination.  Three  pro¬ 
grams  based  on  Gaussian  elimination  have  been  developed. 
These  codes  are  aimed  at  applications  where  the  problem 
is  solved  only  once  and  not  repeatedly  with  the  same 


sparseness  structure.  One  code  uses  column  ordering, 
one  uses  row  ordering,  while  the  third  combines  row- 
column  permutations  to  achieve  reduced  fill-in  and  in¬ 
creased  accuracy.  Also,  a  sparse  matrix  program  tor 
Gauss -Jordan  elimination  was  coded  to  be  used  for  com¬ 
parison  purposes.  . 


Solution  of  Optimal  Control  Problems 


Work  was  continued  on  the  problem  of  the  computation 
of  solutions  of  continuous  optimal  control  problems,  with  the 
goal  of  obtaining  a  procedure  for  solving  problems  with  inter¬ 
mediate  state  constraints.  In  an  earlier  report,  it  was  proved 
that  problems  with  convex  costs,  linear  differential  systems, 
and  any  type  of  control  and  trajectory  constraints  can  be 
approximated  by  discrete  (finite-dimensional)  optimal  con¬ 
trol  problems.  The  proof  given  was  semiconstructive.  These 
theoretical  results  were  tested  on  several  minimum  energy, 
fuel,  and  time  problems  with  intermediate  control  constraints 
and  terminal  state  constraints.  The  discrete  approximations, 
obtained  by  directly  discretizing  these  problems,  were  solved 
by  linear  and  quadratic  programming  algorithms.  The  com¬ 
puted  optimal  costs,  controls,  and  trajectories  were  compared 
with  tlie  solutions  obtained  analytically.  The  results  demon¬ 
strated  that  discretization  is  a  feasible  method  for  these  types 
of  problems. 

One  minimum  time  problem  with  intermediate  state 
constraints  was  also  solved  by  this  technique,  but  the  compu¬ 
tation  was  extremely  inefficient  and  resulted  in  excessively 
large  computational  times.  Effort  is  presently  being  exerted 
to  utilize  the  sparseness  of  the  matrices  in  the  problems 
generated  to  obtain  reasonable  computing  times. 


Time  was  spent  studying  existing  methods  for  tackling 
problems  with  intermediate  state  constraints.  There,  are 
essentially  two  approaches:  (1)  utilizing  the  associated  neces¬ 
sary  conditions  for  optimality;  and  (2)  introducing  penalty 
functions.  Approach  (1)  requires  a  priori  knowledge  of  the 
optimal  trajectory  and  is  very  difficult  even  for  linear  systems. 
Approach  (2)  is  a  direct  method,  but  requires  delicate  treat¬ 
ment  niunerically.  It  was  proved  theoretically  that  approach 
(2),  if  used  to  remove  all  state  constraints,  yields  the  solu¬ 
tion  to  the  problem  to  be  solved  if  and  only  if  the  problem  is 
the  same  as  its  relaxation.  By  example,  it  was  also  demon¬ 
strated  that  if  (2).  is  used  to  remove  only  the  intermediate 


state  constraints,  the  results  obtained  may  be  the  solution  of 
the  original  problem  or  the  solvition  of  its  relaxation.  Hence, 
if  this  technique  is  applied  to  a  problem  that  is  not  the  same 
as  its  relaxation,  care  shoxild  be  exercised.  (Observe  that 
if  the  problem  is  line?  "  with  a  convex  cost,  then  it  is  the  same 
as  its  relaxation.  )  It  is  conjectured  that  this  result  will  be 
true  for  approach  (3)  (direct  discretization),  too. 


Optimal  Design  of  Electrical  Networks 

It  has  been  recognized  for  some  time  that  the 
problem  of  finding  the  optimal  parameter  values  for  a 
given  network  and  a  given  performance  criterion  involves 
the  munerical  calculation  of  the  transient  response  of  the 
network  equations  and  the  adjoint  equations.  With  this 
information,  one  can  compute  the  gradient  direction  of 
the  performance  criterion  with  respect  to  the  parameters. 

The  parameter s, arc  then  changed  in  an  appropriate  way 
and  this  procedure  is  then  repeated.  The  optimum  design 
is  then  obtained  after  many  iterations,  each  one  requiring 
perhaps  1000  to  10,  000  time-step  calculations.  If  the  basic  tim 
step  computation  is  not  efficient,  the  computation  time 
required  can  be  unacceptable.  Since  one  is  involved  with 
stiff  systems  in  electrical  networks  as  well  as  nonlinear 
devices  in  many  practical  design  problems,  the  basic 
time-step  calculation  requires  the  solution  of  a  system 
of  linear  equations,  perhaps  several  times. 

A  preliminary  study  of  this  problem,  starting  with 
the  topological  description  of  the  network,  has  led  to  the 
conclusion  that  network  analysis  can  be  viewed  basically 
as  Gaussian  elimination  for  a  large,  very  sparse,  matrix., 

This  point  of  view  gives  a  uniform  approach  to  network 
analysis  which  allows  one  to  write  a  simple  network 
analysis  package.  The  main  computation  is  based  on 
several  of  the  sparse  matrix  codes  already  developed. 

In  addition,  because  of  its  simplicity,  the  method  allows 
the  inclusion  of  the  most  general  electrical  elements 
(e.  g,  ,  voltage  controlled  voltage  sources)  and  is  easily 
adapted  to  optimal  design  problems.  It  also  frees  one 
of  some  inhibiting  ideas,  about,  how  networks  should  be 
solved,  that  have  developed  in  recent  years.  Since  many 
of  these  ideas  were  developed  without  a  knowledge  of 
sparse  matrices,  they  may  be  inefficient  compared  to  the 
newer  techniques. 

The  conclusions  of  this  study  are  that  a  new  net¬ 
work  analysis  and  optimal  design  package  should  be 


developed,  and  this  is  currently  underway.  In  addition, 
by  modifying  some  of  tlie  sparse  matrix  codes,  to  take 
advantage  of  the  special  nature  of  the  network  problems, 
even  better  efficiencies  can  be  achieved. 
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